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SUMMARY 


Optical-waveguide Luneburg lenses of arsenic trisulfide glass have 
been made and tested. The lenses are formed by thermal evaporation of As 2 S^ 
through suitably placed masks onto the surface of LiNbO^:Ti indif fused wave- 
guides. The lenses are designed for input apertures up to 1 cm and for 
speeds of f/5 or better. They are designed to focus the TMq guided mode of 
a beam of wavelength, external to the guide, of 633 nm. The refractive index 
of the As^S^ films and the changes induced in the refractive index by expo- 
sure to short-wavelength light have been measured. Some correlation between 
film thickness and optical properties was noted. The short-wavelength photo- 
sensitivity has been used to shorten the lens focal length from the as- 
deposited value. Lenses of rectangular shape, as viewed from above the guide, 
as well as conventional circular Luneburg lenses, have been made. Measure- 
ments made on the lenses include thickness profile, general optical quality, 
focal length, quality of focal spot, and effect of ultraviolet irradiation on 
optical properties. The best lenses had speeds of f/3 to f/2, and possessed 
good to very good focal properties — central spots within twice the diffraction^ 
limited intensity width and primary sidelobes no more than 2 db above 
the diffraction-limited intensity — when stopped down to f/10 to f/5. 
Preliminary investigations of information-handling capacity of some lens 
were also carried out. 
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I. INTRODUCTION 


This report describes the results of a one-year research program 
to develop chalcogenide-glass Luneburg lenses for use in LiNbO^-based 
integrated optical devices. More specifically, we describe techniques 
developed in this program, and to some extent in parallel programs, for 
the design, fabrication, and testing of short-focal-length arsenic tri- 
sulfide Luneburg lenses on Ti-indif fused LiNbO^ optical waveguides. As 
part of the program, techniques have been investigated for post-deposition 
modification of the optical properties of the lens by exposure to ultra- 
violet light. Two lenses, one modified as just described, have been 
prepared for delivery to the sponsor as specified in the contract. 

In this report we first outline the theoretical and calculational 
procedures used to design the lenses, to design mask arrangements for 
production of these lenses by evaporation of the glass, and to assess the 
expected lens quality. Following this, we describe the experimental work 
of measuring the required material parameters and of fabricating and testing 
the lenses. In the present section we provide some background Information 
and give a brief summary of our approach in order to motivate and to simplify 
the ensuing discussion. 

A conventional Luneburg lens for use in planar Integrated optics 
consists of a layer of high-refractive-index material deposited on a portion 
of an optical waveguide and radially symmetric about an axis through the 
center of the lens and perpendicular to the waveguide. Figure 1, in which 
this type of lens is illustrated along with the other principal types of 
waveguide lenses, should complement the verbal description. The effective 
refractive index of the guided mode is changed locally by an amount dependent 
on the thickness of the overlay at the point in question. The Luneburg lens 
is an example of a gradient- index lens, in which focusing occurs because of 
the difference in refractive index between adjacent rays. If the lens 
profile, thickness as a function of radius, is properly chosen, a perfect 
geometric focus can be obtained; that is, the lens will be diffraction- 
limited. 
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Fig. 1. Three types of short-focal-length lenses for 
In integrated optics devices 


use 




The lens material should have a refractive index higher than that 
of the waveguide surface index, should not cause excessive absorption or 
scattering of the guided light, and should be easy to deposit on the 
waveguide surface. For lenses on LiNbO^ waveguides, arsenic trisulfide 
(AS2S2) glass is one of the few known materials to meet these criteria. 

complex chalcogenide glasses also have high refractive 
indices, but they have fundamental absorption edges more toward the infrared 
and cannot be used in the visible. ZnS and CdS have large energy gaps and 
high indices of refraction, and they can be deposited on LiNbO^- These 
materials form polycrystalline films; obtaining films of good optical 
quality requires some care. Certain oxides, such as TiO^, can, in principle, 
be used in conjunction with LiNbO^ substrates, but in practice it is 
difficult to keep the oxygen content of the deposited films high enough to 
obtain a refractive index comparable with those of the crystal. 

AS2S2 ^y be deposited either by evaporation from a heated crucible 
or by rf sputtering. This glass has other useful properties as well. It 
may be removed without harm to the waveguide by dissolution in strong 
hydroxide solution, thus minimizing loss of expensive crystal substrates. 

It is also photosensitive, its refractive index increasing substantially 
(by more than 0 . 1 ) on exposure to blue or ultraviolet light. This phenomenon, 
or a similar permanent thermo-optic annealing effect, can be used for post- 
fabrication modification of the lens characteristics. 

A brief survey of the development of Luneburg lenses and of As2S^ 
integrated-optics devices is provided in Appendix A. 

Our present procedure for design and fabrication of Luneburg lenses 
comprises the following steps: 

( 1 ) An ideal Luneburg lens of the desired aperture and focal 
length, is designed, using standard procedures, for the 
waveguide mode and wavelength of interest. Design Involves 
knowing something of the characteristics of the indiffused 
waveguide and selecting a target refractive index for the 
lens material, either the as-deposited index for the design 
wavelength or that expected after post-deposition modification. 
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(2) A mask arrangement suitable for producing a good 
approximation to the design lens profile is devised. 

In our work so far the lenses have been made by evaporation 
of the glass from a crucible through the mask apertures 
onto the waveguide; so the size of the crucible must be 
considered in the design. We have designed principally 
with only two masks, one to obtain the main part of 
the lens curvature and one to adjust the size and contour 
near the edge. The best mask positions and apertures 
are found approximately by a coarse optimization 
procedure and then refined by trial-and-error calculations. 

(3) The lens is formed, as mentioned above, by evaporating 

powdered ^^ 2^3 from a quartz crucible onto the 

waveguide. The evaporation is carried out in a vacuum 
of around 1.0 x 10*^ torr. A piezoelectric thickness 
gauge is used to measure the amount of material deposited. 
The evaporation and masking configuration is sketched in 
Figure 2. 

(4) If the lens is to undergo post-fabrication modification, 
it is subjected to a controlled exposure to ultraviolet 
light; otherwise it is protected from ultraviolet. 

After a lens is made, it may be tested, optically and physically, 
in a variety of ways which we will describe in the appropriate place in 
the report. 


II. THEORETICAL AND CALCULATIONAL WORK 
A. Lens Design 

The first step in lens design is, of course, selection of desired 
lens characteristics. Most of the lenses discussed in the present report 
have been designed for a focal length, within the waveguide, of 30 mm and 
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Lens (Height 



Source 


^^8* 2. Model of experimental arrangement for evaporation of AS 2 S 2 Luneburg 

lenses. Source is actually a crucible filled with powdered As-S- glass 



an input aperture of 10 nnn. They are intended to focus a TM^ mode, of 
effective refractive index 2.288 in the waveguide region external to the 
lens. The designs are for light of free-space wavelength 633 nm. This 
set of characteristics provides a rather stringent test of design and 
fabrication capabilities. 

With just the information given so far, we can determine the 
refractive index profile required for the lens by solving the Luneburg- 
Morgan integral equation: 


In 


N(r) _ 2 


(t - z) 


1/2 


•w 


ext 


sin ^ (u^ + z) 


du. 


( 1 ) 


In this equation, N(r) is the mode index at radius r in the lens region, 
while is the mode index outside the lens. The parameter t is the 

reciprocal of twice the f /number of the lens, and z = tR, where 
R = 2r N(r)/S Next. A is the lens aperture. In the form presented here, the 
integral is easily evaluated to 5 decimal place accuracy by a single 
16-point Gaussian quadrature. The design refractive index profile is 
shoxm in Figure 3. 

To find the lens thickness profile corresponding to this refractive- 
index profile, we need the refractive-index profile of the diffused wave- 
guide and the bulk refractive index of the lens material. We have taken 
the waveguide to have a Gaussian profile with a diffusion length of 2 ym, 
a surface index of 2.293, and a bulk index of 2.287. These values, which 
are not highly critical to the accuracy of the calculations, are typical 

of single-mode LiNbO^rTi waveguides for the polarization and wavelength 
conditions stated. Since the lens thickness varies slowly with radius, 
we can model the situation at a given radius as a uniform layer of lens 
material covering the inhomogeneous waveguide. We made a straightforward 
extension of the calculation method devised by Southwell (1977)* to this 


* See the bibliography in Appendix A for citations of work referred to 
by author and year. 
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case. This analysis applies, strictly speaking, only to TE modes, except 
of course that index values appropriate to the TM^ mode in the LiNbO^ are 
inserted in the calculation. Subsequent work [see also Bryan et al. (1982)] 
has shown that while this is in most circumstances a fairly good approxima- 
tion, it should be corrected for more accurate lens design and to allow 
better identification of sources of observed aberration; the required simple 
program modifications have not yet been made. 

The refractive index at 633 nm of the AS 2 S 2 material was measured 
by methods described elsewhere in this report. We found the index of as- 
deposited ^^ 2^3 about 2.445, with a fair amount of variation 

from sample to sample, which we presently believe may result mainly from 
variations in film stoichiometry with thickness. We made lens designs for 
film indices of 2.445, 2.46, and 2.48. The relative lens profiles, normalized 
to unity at the center of the lens, are shown in Fig. 4. The actual design 
central thicknesses are 1.69, 0.83, and 0.55 ym, respectively. This marked 
variation of required film thickness with assumed film index, and conversely 
of required film index with actual film thickness, brings out the importance 
of the post-fabrication index-modification capability. 

B. Mask Design 

Given the lens profile required for a particular focusing require- 
ment, we need to determine apertures and positions of masks suitable for 
depositing a lens with this profile on the waveguide. We have found that 
fairly good approximations to profiles of interest can be obtained using 
only two thin masks. 

The model for the evaporation system, and the system itself, have 
already been illustrated schematically in Figure 2. 

The evaporation source is modeled as a uniform distribution over 
a circle of radius b of point sources, each emitting AS 2 S 2 molecules uniformly 
into the hemisphere above the source. The extent to which this simple 
model describes actual evaporation from a crucible has not been investigated 
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relative thickness 



lens radius, mm 


Fig. 4 . Design profiles for As2^2 huneburg lens on LiNbO^ waveguide supporting a 
single TM mode for three values of lens refractive index n^. 

(a) rif - 2.445; (b) n^ = 2.46; (c) n^ = 2.48. 
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experimentally. It is possible that the geometric radius b' of the top of 
the crucible or the top of the melt differs from the effective source radius 
b that must be used in the model. 

The lens-profiling mask is modeled as a simple circular hole in 
an infinitesimally thin sheet parallel to the waveguide surface. The 
distance h^^ of this sheet from the surface and the radius s of the hole 
are varied, first by a coarse least-squares procedure at five or six points 
and then by trial and error to optimize the fit to the desired profile. 

The fits are to a lens profile normalized to unity at the center; we rely 
on a separate measurement of the central thickness to obtain the proper 
lens thickness. To obtain a good fit near the edge of the lens, an additional 
mask very close to the substrate is used. It is positioned to reduce the 
lens thickness to zero at the desired aperture. 

One cannot obtain lens profiles with regions of substantial upward 
curvature, such as shown in curve (a) of Figure 4, with only two masks. We 
tried nonetheless to make with two masks as good a fit as possible to this 
profile in order to test the adequacy of this simple arrangement. The 
expected attainable profile is compared to the design profile in Figure 5. 

The mask conditions to obtain this profile are as follows (dimensions in mm) : 

Mask Distance below sample Diameter 


Source 

100 

18.3 

Profiling 

28 

5.0 

Edge-def ining 

0.5 

11.9 


Although the lens profile appears rather steep at the edge, lenses 
made according to this prescription did not produce any noticeable scattering 
into higher-order modes. Such scattering is readily observed since it 
produces overlapping output beams with two or more different focal lengths 
and consequently no single sharp focal plane. It is Important that the 
overlay film be sufficiently tapered in thickness at the lens edges that 
these spurious modes in the lens region not be excited. Such modes occur 
when the required overlayer thickness is sufficient to allow higher modes 
in the underlying waveguide to exist. As the film thickness is increased. 
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the single mode of the waveguide is slowly "pulled" up in effective index, 
until a point is reached when the index rapidly increases beyond the wave- 
guide surface value, toward the bulk index of the film. As this rapid change 
is occurring, a new mode appears in the waveguide, a mode that was cut off 
for thinner films. This mode replaces the original single mode of the 
guide; its index rapidly increases to the vicinity of that of the original 
guided mode, where it remains until it too is finally pulled up into the 
film. This process continues, with new modes being pulled up into the 
guide from cutoff to replace lower modes that have been pulled into the 
film. As a result, there is, roughly speaking, almost always a mode present 
having an effective index very near that of the original waveguide. If 
the film thickness changes too rapidly at the edge, then in the region 
where the original mode is lifted into the film, strong coupling can occur 
between the incoming mode and the spurious mode lifted from the cutoff by 
the presence of the film. This coupling occurs because the effective index 
difference between the incoming mode and the spurious one is small and is 
changing more slowly than that between the incoming mode and its continuation 
in the film region. The result of such coupling is the appearance of 
multiple focused spots and increased apparent loss of transmission. The 
increased loss may not be too serious, but the transmission of the spurious 
modes through the lens is deleterious for most applications and can be 
disastrous for optical data processing. 
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III. EXPERIMENTAL WORK 


A. Parameter Measurements 

For accurate lens design and modification, and for calculation of 
lens properties, it is necessary to know several physical parameters of the 
lens material; consequently we have devoted considerable effort to measuring 
these parameters as determined under our evaporation conditions. 

1. Lens Shape . As mentioned previously, the overall thickness of 
the lens is fixed by measuring the mass of material deposited on a piezo- 
electric thickness monitor mounted near the crystal but away from the masks. 
The mass/thickness ratio is determined by profilometer thickness measurements 
on test specimens. The profilometer, a Talysurf 4 manufactured by Taylor- 
Hobson, is also the principal Instrument used to determine the lens profile. 
The stylus was not found to damage the arsenic trisulfide glass or the 
waveguide layer, and it is not too difficult to make a traverse through the 
thickest part of the lens. One difficulty with the Talysurf measurements 
is that the substrate near the lens was often found not to be flat; so it 
was necessary to make an extrapolation of substrate surface position beneath 
the lens in order to determine the lens thickness. It seems that the 
deposition of the lens may cause or enhance some warping of the substrate, 
especially as the effect seems to be more marked on thinner LiNb 02 crystals. 
Departure of the substrate from planarity was not, however, observed to be 
sufficient to have any effect on the properties of the guided waves or of 
the lenses. 

The profilometer work is complemented by interferometry. A Twyman- 
Green arrangement is used to provide interferograms of the lens shape; an 
example is shown in Figure 6. Both white-light and monochromatic illumination 
yield informative fringe patterns. Quantitative analysis of such fringe - 
patterns has been limited to confirming approximate central thickness of 
the lens by fringe counting; but what the interferograms are particularly 
good for 'is detecting shape distortions resulting from misalignment of 
masks or from substrate imperfections. 
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Fig. 6. Twyman-Green interferogram, made with 488 nm light ^ 
of As 2S-:^ Luneburg lens, 1.2 ym thick at center and 
1 cm in diameter, on LlNbO^ substrate. 



The refractive index at 633 nm of 


2. As«S^ refractive index. 

— I j 

arsenic trisulfide films deposited according to our procedure has been 
measured by a prism deflection technique. Prisms of uniform thickness and 
apex angle 30® were deposited on LiNbO^ waveguides through triangular masks 
held close to the guide. The prisms were oriented symmetrically with respect 
to the z-axis of the crystal and a 633 nm guided beam was coupled into the 
waveguide so it propagated along the x-axis of the y-cut crystal. After 
the beam was deflected by the overlay prism, it was end-fired out of the 
guide and its deflection was measured on a screen about a meter away. 

Because of the steep edges of the prisms, some of the light in the prism 
region was frequently scattered into, and subsequently back out of, higher- 
order modes, so there were usually 2 or 3 deflected beams corresponding to 
different modes in the prism region. For a 30® prism, the x-directed input 
beam exits at close to the expected angle of minimum deviation; it is 
difficult to find and to analyze the minimum deviation conditions accurately 
though, because of the variety of media through which the beam passes; so 
we analyzed the simpler situation of the initially x-propagating beam. From 
the observed deflection we calculate the mode indices in the prism region 
and from each mode index we calculate a value for the refractive index of 
the as-deposited overlay material. The unweighted average of 16 such 
determinations, measured on 6 prisms varying between 0.28 and 1.21 ym in 
thickness, was 2.445 + 0.006. Two different substrates were used, and the 
modes observed were 6 TEq, 6 TMq, 3 TE^ and 1 TMj^. The range of refractive 
indices found was 2.38 to 2.48, larger than desirable, but most of the 
values clustered well around the mean; consequently we adopt 2.445 as a 
provisional value in the design work. For the 7 TM modes the mean was also 
2.445, while for the 6 TM^ modes it was 2.447. 

3. Ultraviolet-induced change in refractive index . The dynamic 
photoinduced index change was determined by measuring the change in the 
deflection of a guided TM^ beam by an As^S^ prism as the prism was illinninated 
with an ultraviolet lamp. This prism again had apex angle 30®; it was 
0.60 ym thick. The single-mode waveguide on which it was placed had polished 
end faces to improve the quality of the beam end-fired from the guide. The 
input coupling prism and the crystal were all enclosed in a plastic box 
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through which dry argon was flowed throughout the experiments. The ultra- 
violet illumination, strongest at 400 nm, impinged on the prism through the 
top of the box; the intensity of the 400 nm line at the sample, with the 
box cover in place, was 1.10 mW/cm^. It was noticed that there was a small 
change in refraction by the prism between the time it was first made and 
when it was used in these experiments, although the sample had been kept 
in a desiccator and under yellow light. Backward extrapolation of the 
photoresponse dynamic curves indicated that this preliminary change in film 
index was equivalent to about 10 minutes extra exposure to the ultraviolet 
source. For exposure times up to about 3 hours, it was fovind that the 
refractive index of the ^ 32 ^^ as measured by the deflection of the TMq mode 
could be expressed as 

n = 2.5820 - 0.1295 exp (-E/5.739), (2) 

2 

where E is the total uv exposure (400 nm line) in J/cm . Such an expression 
is consistent with a model in which the rate of photoinduced index change is 
simply proportional to the remaining amount of unaltered material. For 
longer exposures up to 8 hours the film index increased above the saturation 
value indicated by Eq. (2) to 2.594. This effect is probably related to 
heating of the sample by the uv lamp; when the refractive-index changes 
were becoming small, long exposures (1 hour) were made between measurements. 
The prism region also supported a TM^^ mode; the dynamics of the index change 
measured using deflection of this mode were very similar to those with the 
TMq mode, but the apparent initial film index calculated for this mode was 
about 0.02 higher than that found using TMq. After 6 hours exposure a TM 2 
mode of even higher apparent index showed up. The observed index changes 
and model curves fitted to an equation of the type of (2) are shown in 
Figure 7. While heating of the sample seems to have a major influence on 
the values obtained at long times, the difference in the Initial values 
should not go unremarked. We presently believe this difference may result 
from some dependence of film composition on evaporation time. Tanaka (1980) 
reports a fairly sharp peak in the ultraviolet photosensitivity at the 
stoichiometric composition The refractive index of the as-deposited 
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0 100 200 300 400 500 


ultraviolet exposure time, minutes 

Figure 7. Apparent refractive index of AS 2 S 2 film prism as function of time of exposure to 1.1 mW/cm^ of ultraviolet 
light, as determined by deflection of 3 different modes supported in prism region. Heavy lines connect 
data points; light lines are fits to data of the form n^ = a - b exp(-t/x). 




film is also maximtim near this composition; it falls off at higher sulfur 
concentrations but decreases little if at all, at more arsenic-rich 
compositions. It seems probable, therefore, that our films are near- 
stoichiometric when thin, but become somewhat sulfur-rich when thick. The 
vapor is believed to be sulfur-rich at all times during the evaporation; so 
the differences in film composition may be related to differences in 
condensation properties between the lithium niobate substrate and the 
chalcogenide glass itself. Higher-order guided modes, which have a mode 
index closer to that of the substrate than does the fundamental, yield higher 
film indices calculated on the basis of a uniform film because these modes 
are more effectively drawn into the high- index region near the surface. 
Several other experimental results, described elsewhere, tend to support 
this interpretation. 


bens Fabrication and Testing 


In this section we provide a few details regarding our present 
waveguide and lens fabrication and testing procedures. 

o 

► The optical waveguides were formed by evaporating 175 A of 
titanium onto the surface of Y-cut LiNbO^ crystals 1 to 3 mm thick and 
indif fusing this film at 1000 C for 2.5 h in an atmosphere of flowing 
oxygen which had been bubbled through water at 90 C. This procedure 
produces good quality waveguides supporting a single mode of each polariza- 
tion at the design wavelength of 633 nm. 

The lenses are made in a conventional bell- jar high-vacuum 
evaporation system, typically pumped down to a pressure of about 1.0 x 10 ^ 
torr at the start of the evaporation. The source is 99.99% pure AS 2 S 2 
fused glass which has been hand-ground in a porcelain mortar and pestle 
to a fine powder, approximately 325 mesh. The powder is evaporated from a 
quartz crucible held in a tungsten basket. The evaporation temperature is 
estimated to be in the 500 to 700 C range and the source-to-substrate 
distance is typically 100 mm. In our earlier experiments, the film 
deposition rate was around 20 nm/second. Since we considered that this 
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rapid evaporation rate might be a source of change in lens composition with 
thickness, we have more recently fabricated lenses using deposition rates 
of 1 to 2 nm/sec. While these more slowly deposited lenses are generally 
of better quality, there is so far no evidence that they are of more 
homogeneous composition. The masks are made of thin sheets of aluminum. 
The crystal substrate was not heated or cooled during the evaporation. 

Post-fabrication modification of the lens refractive index was 

generally carried out with a more intense uv source than that used in the 

study of the dynamics of the index change. The intensity at the sample, 

2 

measured at the dominant 400 nm line, was usually around 20 mW/cm . 
Reciprocity has been assumed in estimating index changes resulting from 
exposures at this level. 

Our primary method for characterizing the optical quality of the 
Luneburg lenses has been the examination of the light distribution in the 
focal plane. The focal spot was scanned by coupling the beam transmitted 
through the lens out through a rutile prism and refocusing it with an f/2 
imaging lens onto an optical multichannel analyzer. The experimental 
arrangement is shown in Figure 8. The OMA has 500 25 ym channels on 25 ym 
centers. The channels are long enough to collect substantially all the 
light diffracted or scattered in the direction perpendicular to the wave- 
guide plane. In earlier work, we used a diode array or a fiber scanner to 
investigate the focal spot; in these cases an additional cylindrical lens 
was necessary to collect all the light. The OMA output can be displayed 
on an oscilloscope screen or recorded digitally. 

In most cases the lenses focused outside the waveguide. To 
determine the focal length, we measured the length of the optical path in 
each medium — waveguide, output coupling prism if used, and air — from the 
lens center to the focal point, and converted the total distance to an 
equivalent distance inside the waveguide. The accuracy of such a procedure 
clearly cannot be very high; and our quoted focal lengths might easily be 
in error by +10%. 

The information-handling capacity of the lenses has so far been 
investigated in only a preliminary way by placing a Ronchi ruling in the 
input beam and observing the relative intensity and width of the sidelobes 
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Optical Multichannel 

Variable Analyzer or Diode 



Fig. 8. Experimental arrangement for determining quality of lens focal spot. 

When the optical multichannel analyzer is used, the cylindrical lens 
is not needed. 
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generated in the output pattern. Only a few different rulings have been 
available. A better procedure seems to be to place some lenses on wave- 
guides also containing broad-band surface-acoustic-wave transducers. These 
can be used to generate phase gratings of a range of spatial frequencies in 
the lens input plane, which should permit a determination of minimum and 
maximum resolvable spatial frequencies of the lens. 

C. Summary of Some Experiments 

Numerous lenses have been fabricated and tested. Some lenses 
showed a high degree of scattering of the guided light, so were not 
examined further. Usually the scattering appears to be related to the 
presence of trace amounts of diffusion-pump oil; thus, this source of 
scattering can presumably be eliminated by more frequent system maintenance. 
Some design, physical, and optical characteristics of several lenses are 
given in Table I. It will be noticed that the lens central thickness is 
frequently markedly different from the design thickness. In some cases 
this results from failure of the thickness monitor; in others it is 
intentional, to investigate the effect of varying lens thickness on the 
optical properties. The initial focal length does not vary as expected 
with the lens thickness. This may reflect changes in stoichiometry with 
lens thickness as previously discussed. 

OMA scans of the lens focal spots have generally been made at 
reduced input aperture in order to be sure that all the light transmitted 
is captured by the relay lens and focused on the OMA detector. The 
diffraction patterns do not sharpen markedly at larger apertures, but we 
are at present not able to say how much of this effect, or noneffect, 
results from poorer quality of the lens near the periphery and how much 
from aperture effects in the light-collection system. The focal spot 
quality does not appear to vary in any marked or predictable way with lens 
thickness. The spot quality data may be compared with the ideal by using 
Eq. (1) of Appendix B. For 4 mm aperture and 25.6 mm focal length, the 
half power diffraction-limited spot size for the TMq mode at 633 nm is 1.6 ym. 
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TABLE I. PROPERTIES OF AS-PREPARED LENSES 


Design 


Measured 


N) 

u> 


Sample 

Aperture, mm 

Mask 

Thickness, 

Focal Length, 

Spot Size, 

Sldelobe, 
height, dB 


Number 

useful/total 

Arrangement 

pm 

mm 

pm 

Comments 

200A 

6/8 

(a) 

0.36 

42.8 

8.32 

-9.4 

- 

198A 

6/8 

(a) 

0.30 

168.7 

- 

- 

- 

199 

6/8 

(a) 

- 

56.5 

- 

- 


200 B 

6/8 

(a) 

- 

125.2 

- 

- 

Samples 200 & 198 








were reused. 

198 B 

6/8 

(a) 

- 

273.0 

- 

- 

- 

124 

6/8 

(a) 

- 

43.0 

- 

- 


HP-6 

10/12 

(b) 

0.68 

IS.O(TEo) 



Subsequently annealed 





28.7(TM^) 

2.9 

-7 ^ 

at 190 C for 30mln. 
in flowing Argon. 
Focal length re- 
duced to 15.0 mm 
(TMo) 

250 

10/12 

(b) 

0.96 

22.2 

_ 

_ 

As deposited. 

166 

10/12 

(b) 

(b) 

1.20 

26.7 

- 

- 

thicker lenses have 
longer focal lengths 

226A 

10/12 

1.20 

31.3 

- 

- 


206 

10/12 

(b) 

1.46 

31.3 

- 

- 

- 

253 

10/12 

(b) 

0.80 

17.6(TEo) 

5.8 

-11.1 




(b) 


22.2(TMjj) 

5.5 

-6.2 


226 B 

10/12 

(b) 

0.70 

17.6(TEq) 



Sample #226 was re- 
used. TE^ focal 





24.5(TMq) 

2.28 

-11.2 

spot was also sharpy 








1 but was not scanned 

W-1 

10/13 

(c) 

0.56 

92.3 

7.3 

-4.3 

Rectangular 







(4x13 mm) 

W-2 

10/13 

(c) 

0.89 

115.2 

8.84 

-4.8 

Rectangular 


*(a) Source to LlNbO^ 100 mm; profiling mask aperture 4.1 mm, position 23.5 mm below waveguide; 
edge-defining mask aperture 7.8 mm, position 0.75 mm below waveguide. 

(b) For mask arrangement see p. 11. 

(c) For mask arrangement see p. 44. 



For the ideal lens the first sidelobes should be 13.3 dB down in intensity 
from the central peak. The two numbers given, focal spot size and first 
sidelobe intensity, do not fully characterize the lens quality, of course, 
even at fixed input aperture. It should also be borne in mind that the 
input beam is a truncated Gaussian rather than an ideal plane wave. In 
Table II we summarize effects of short wavelength light, from a variety 
of sources, on the optical properties of a number of lenses. While 
ultraviolet exposure reduces the focal length markedly by increasing the 
refractive index of the lens material, it does not generally lead to 
marked deterioration in apparent focal spot quality. After long exposure 
the fraction of light in the sidelobes tends to increase, though; and 
sometimes, as in lens W-2, it becomes quite large. 

Further lens testing was done by placing Ronchi rulings in the 
input beam and examining the resulting focal pattern. The experimental 
setup is sketched in Figure 9. The light passing through the lens was 
.coupled out of the waveguide and the focal plane intensity imaged at a 
magnification of 44. Photographs of such re- imaged focal planes for 
lens 226B are shown in Figure 10. The principal diffraction peaks from a 
grating of 2.36 lines/ram (60 lines/inch) are imaged adequately, and those 
from a grating of 7.87 lines/mm (200 lines/inch) are well-separated. The 
positions of these peaks were measured on* the original photographs; they 
occurred where expected. For an ideal Ronchi ruling illiminated by a 
plane wave, the second-order principal diffraction peak should be absent. 
Traces of this peak are apparent in both photographs, although with the 
coarser ruling it is evident that the third-order peaks are stronger. 
Further experiments to characterize more fully the spatial information 
handling capacity of the lenses are planned. 



TABLE II. EFFECTS OF EXPOSURE ON LENS PROPERTIES 


Sample 


Focal Length, 

Spotslze, 

- At Aperture, 

Side lobe 


Number 

Exposure- Source 

mm 

\sa 

mm 

Helght,dB 

Comments 

200A 

0 

42.8 

8.32 

4 

-9.4 

.. 


0.25J-25W bulb^^^ 

24.5 

4.87 

4 

-11.1 

- 

198A 

0 

168.7 

- 

- 

- 

- 


0.5J-25W bulb 

26.7 

6.32 

4 

-13.0 

- 

199 

0 

56.5 


_ 

1 

Exposure measured with 


0.60J-Argon Laser 

45.1 

- 

- 

i 

broadband laser power 
meter 


2.4J-Argon Laser 

38.2 

- 

- 

1 


200 B 

0 

125.2 






eOmJ-Blackllght^*^^ 

104.6 

- 

- 

- 

Exposure measured with 
broadband laser power 


90mJ-Blacklight 

93.2 




meter 


150mJ-Blackllght 

81,7 






540mJ-Blackllght 

54.2 





198 

0 

273,0 

- 

- 

- 



0,13J-25W Bulb 

103.0 






0.20J-25W Bulb 

63.0 






0.46J-25W Bulb 

43.0 





124 

0 

43.0 






2. 16 J- Flood Lamp^^^ 

38.0 






14 • 4 J- Flood Lamp 

33.0 







K) 


TABLE II. (Continued) 


Sample 

Number 

Exposure-Source 

Focal Length, 
mm 

Spot Size - 
pm 

At Aperture, 
mm 

Sldelobe 
Height, dB 

Comments 

250 

0 

22.2 

_ 

_ 

- 

- 


61J-UV Lamp^®^ 

17.6 

- 

- 

- 

- 


122J-UV Lamp 

13.0 

- 

- 

- 

- 

166 

0 

26.7 

- 

- 

- 

- 


61J-UV Lamp 

22.2 

- 

- 

- 

- 


122J-UV Lamp 

17.6 (TM^) 
17.6 (TE^) 

2.93 ym 

4(TMo) 

-13.6(TMo) 

Very good quality lens 

226a 

0 

31.3 

- 

- 

- 

- 

206 

61J -UV Lamp 
122J-UV Lamp 

0 

61J-UV Lamp 

22.2 

17.6 

31.3 

26.7 

- 

- 

_ { 

Both samples 226A and 
206 did not focus very 
sharply as deposited, 
but did after UV 
exposure. 


122J-UV Lamp 

22.2 





W-1 

0 

92.3 

7.3 

3 

-4.3 

Rectangular 


1.98J-UV Lamp 

92.3 

11.34 

3 

-6.3 

Spot size increased, side 
lobes decreased 


32.6J-UV Lamp 

56.6 

5.35 

3 

-5.1 

Optimum exposure 


93.8J-UV Lamp 

51.0 

5.91 

3 

-3.9 

- 

W-2 

0 

115.2 

8.84 

3 

-4.8 

Rectangular 


30.6-UV Lamp 

92,3 

9.18 

3 

-5.3 

Spotsize increased, side 
lobes decreased 


91.8J-UV Lamp 

80.8 

7.18 

3 

-3.1 

Spotsize decreased, but 
side lobes increased 


153.0J-UV Lamp 

69.4 

6.95 

3 

-2,7 




TABLE II. (Continued) 


25W incadescent bulb, measured at 400 nm, + 50 nm bandwidth. 

488, nm line of Ar-lon laser 
Cc ) 

Low Intensity long-wavelength-UV blacklight. 

125W floodlamp, measured at 400 nm. 

(e) 

Ultraviolet lamp, measured at 400 nm. 


to 




Variable 



Fig, 9. Experimental arrangement for investigation of spatial-frequency 
response of lenses using Ronchi rulings. The input beam must be 
well collimated. 


to 

00 



no ruling 


2.36 lines/mm 
Ronchi ruling 



7.87 lines/mm 
Ronchi ruling 


Focal plane images of Ronchi rulings placed 
in input beam to lens 226B 





IV. OTHER ACTIVITIES 


A. Rectangular Lens 


Integrated-optics Luneburg lenses as conventionally described have 
a circular outline as viewed from above the waveguide. Other shapes are not 
prohibited, though; in particular, W. H. Southwell had described Luneburg- 
type lenses of rectangular outline which promise considerable savings in space 
on the waveguide when a very large field of view is not required. We have 
made, by evaporating AS 2 S 2 onto LiNbO^;!! waveguides, some lenses of this 
type and measured their optical characteristics. Our work is described in a 
paper to be submitted for publication shortly. A preliminary version of this 
paper is included as Appendix B. 

B. Presentations and Publications 

Talks describing some of the work on this program were given by 
J. R. Busch at the NASA Conference on Optical Information Processing for Aero- 
space Applications, August 1981, and by V. E. Wood at the Optical Society of 
America Annual Meeting, October 1981. Information on publications resulting 
from these talks, and on the planned rectangular lens paper, is summarized 
in Appendix C. 
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V. CONCLUSIONS 


This program has been successful in showing that short-focal- 
length lenses of good quality can be reliably fabricated by depositing 
arsenic trisulfide through suitable masks onto LiNbO^sTi waveguides. A num- 
ber of incompletely resolved questions remain, however, and several process 
improvements will be required if near diffraction-limited performance is to 
be reproducibly achieved. Among the areas where further work is desirable 
are: improved control of composition and homogeneity of the lens material, 

minor improvements to lens design programs, more general mask design, ray- 
tracing studies of fabrication tolerances, and refinement of photosensitivity 
data. 
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APPENDIX A 


HISTORICAL SKETCH 




APPENDIX A 


HISTORICAL SKETCH 


Lenses of the type associated with the name of R. K. Luneburg were first 
described by him in lectures of 1944. His discussion covers more general 
lenses than the sphere with n2(r) = 2 - r2 mentioned by Born and Wolf. 
Generalizations of Luneburg's ideas were made by a number of workers; papers 
by Kunz (1954), Toraldo di Francla (1955), and Morgan (1958, 1959) summarize 
much of this early work, which was mostly directed toward microwave 
applications. The Luneburg lens was introduced into integrated optics 
by Zernike (1974), who deposited Nb 205 lenses through a tapered mask onto 
glass waveguides. Design principles for integrated-optics Luneburg lenses on 
uniform waveguide layers were given by Southwell (1977a, b). 

Yao and Anderson (1978, 1979) demonstrated diffraction-limited f/5 Luneburg 
lenses made by sputtering Ta 205 onto uniform glass waveguides. 

Hatakoshi and co-workers (1979) described f/2 lenses formed by evaporation of 
ZnS onto glass waveguides. Both Yao et al and Hatakoshi et al discuss briefly 
the design of profiling masks. Their lenses were designed to focus a TE mode, 
to which Southwell's analysis applied. Bryan et al (1982) adapted this work to 
the design and fabrication of TM-mode Luneburg lenses of Ta 205 on glass 
waveguides. Luneburg lenses on inhomogeneous waveguides, such as are formed by 
diff fusion or by ion exchange, were first described by Busch et al (1981), 
who also discussed the modifications to Southwell's design procedure necessary 
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to treat this situation. Recent topics of interest include development of more 
tolerant designs for given depth of focus ( Colombini , 1981); the possibility 
of ion-implanted Luneburg lenses ( Sochacki > 1982); and lenses deviating from 
circular symmetry ( Southwell, 1981, 1982). 

The development of AS 2 S 3 glass has been reviewed by Frerichs (1953). 

^^ 2^3 films were first considered for use in integrated optics because of 
their good acoustooptic properties ( Dixon, 1967; Ohmachi , 1973), although they 
display high acoustic attenuation. Young (1971) was among the first to prepare 
evaporated glassy AS 2 S 3 films of good optical quality. The use of such 
films for waveguiding and for other integrated-optics applications was promoted 
by Klein (1974), The photo-induced changes in refractive index were described 
by Keneman (1971) and by Ohmachi and Igo (1972). Related changes in optical 
absorption were first noted by Kostyshin et al (1966) and were described more 
completely by Berkes et al (1971) and by Igo and Toyoshima (1972). Chemical 
changes, revealed by etching-rate experiments, related to the optical changes 
were found by Keneman (1974). Films made by rf sputtering were first reported 
by Watts et al (1974). The dependence of the photo-induced changes on glass 
composition has been extensively studied by Tanaka and co-workers; a recent 
paper (Tanaka, 1980) contains references to earlier work. A number of models 
have been proposed to explain the photoeffects at an atomic or molecular level; 
papers by Halpern (1976), Bishop and Taylor (1979), 

and Biegelsen and Street (1980) provide physical descriptions of several of 
these models. Among the many integrated-optical devices made using AS 2 S 3 
films, we mention only a multiplexer ( Handa et al, 1980) and a 
light-controlled switch ( Tanaka and Odaj ima, 1981). The remarkable 
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glass-forming tendency of AS 2 S 3 is discussed by Phillips (1982), 

Many additional works might be cited; this brief survey has dealt 
principally with those bearing on the present study. We have also included a 
number of references which survey other work and provide further references, A 
bibliography of works we have referred to follows. 
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RECTANGULAR LUNEBURG-TYPE LENSES FOR INTEGRATED OPTICS 

Van E, Wood and James R. Busch 
Battelle Columbus Laboratories, 505 King Avenue, Columbus, 

Ohio 43201 

Duncan T. Moore and Ben Woolley 
Institute of Optics, University of Rochester, Rochester, 

New York 14627 

W, H. Southwell 

Rocketdyne Division, Rockwell International, Canoga Park, 

California 91304 


Abstract 


Compact Luneburg-type lenses of rectangular outline as viewed from 
above have been made by thermal evaporation of AS 2 S 3 glass onto 
single-mode LiNb 03 :Ti waveguides through suitably shaped masks and 
subsequent exposure of the glass to ultraviolet light. The best lenses 
had speeds of f/5,5 at an aperture of 10 mm and focal spots at reduced 
aperture about 1,2 times the diffraction-limited size. These lenses have 
a field of view of at least 25 degrees inside the waveguide. 


Among the types of lenses which may be used in integrated optics 
devices, Luneburg lenses possess advantages of ease of fabrication and 
freedom from aberrations. Diffraction-limited Luneburg lenses have been 
made (0 by sputtering Ta 205 onto glass waveguides; lenses of fairly 
good quality made by evaporation of AS 2 S 3 glass onto Ti-indif fused 
LiNb 03 waveguides have also been described v 2)^ Conventional Luneburg 
lenses (3, 4) ^^e circular in outline as viewed from above the waveguide 
and thus have the additional advantage of a 360 degree field of view. If 
an unlimited field of view is not required, lenses of different outline 
are conceivable; in particular, it has recently been shown(5) that 
diffraction-limited Luneburg-type lenses with a rectangular outline can be 
designed. These lenses offer considerable savings in space, particularly 
important on single-crystal substrates. They may also offer advantages 
over conventional Luneburg lenses of wider fabrication tolerances and 
reduced probablility of scattering into the higher-order guided modes 
which are usually present in the lens region. In the present paper, we 
describe Luneburg-type lenses of approximately rectangular outline formed 



by theraal evaporation of AS2S3 through suitable masks onto single-mode 
LiNb03 waveguides. 

O 

The optical waveguides were fomed by evaporating 175 A of titanium 
onto the surface of a Y-cut LiNb03 crystal 3 mm thick and indiffusing 
this film at 1000 C for 2.5 h in an atmosphere of flowing oxygen which had 
been bubbled through water at 90 C. This procedure produces good quality 
waveguides supporting a single mode of each polarization at a wavelength 
(in air) of 633 nm; Nomarski micrographs show a surface free of ‘’orange 
peel” or other indications of residual oxides of titanium or other metals. 

No formal design procedures have been followed so far in making the 
lenses ; rather we have drawn on experience with conventional 
circular-outline Luneburg lenses and also have relied on trial and error; 
unsatisfactory lenses are easily removed without damage to the waveguide 
by dissolving them in KOH. The lenses are formed by evaporating powdered 
AS2S3 glass from a crucible of circular cross-section, 18.3 mm in 
diameter at its mouth and 110 mm from the waveguide surface. The crystal 
substrate was not heated during the evaporation. The lenses were 
configured using one circular-aperture and one rectangiilar-aperture mask. 
For example, one lens was made by placing the circular mask, which had an 
aperture of diameter 5.0 mm, 38.0 mm below the crystal; and the 
rectangular mask, which had an aperture 2 nm by 11 mm, 10.0 mm from the 
crystal. After the deposition, the lens was Illuminated with an 
ultraviolet lamp to increase its refractive index. The total exposure, 
measured at the dominant 400 nm line of the lamp, was 33 J/cm2. 

According to prism-deflection experiments we have carried out, this 
exposure should increase the glass refractive index to about 2.58. A 
Twyman-Green interferogram of this lens is shown in Fig. 1. The large 
fringes visible on the periphery of the lens arise from flexing or bowing 
of the LiNb03 substrate. The lens Is^ 13 mm wide and 4 mm thick. 

Talysurfs of the two principal directions showed that it is about 560 nm 
high at the center; this measurement was confirmed by fringe counting on 
the interferogram. 

The focusing properties of a similar lens for a 9 mm wide truncated 
Gaussian input beam propagating in the TEq mode are shown in Fig. 2. The 
light source was the 633 nm line of a He-Ne laser. No effects of 
scattering into higher-order guided modes in the lens region were noted. 
The measured focal length is 35 mm for the normal-incidence beam; so the 
lens speed is f/3.9 for the aperture used or f/2.9 for the full lens width 
(12 mm). The useful field of view of the lens, as judged by the range of 
input directions which yielded an apparent sharp focus for an input beam 9 
mm wide, is at least 25 degrees, as measured inside the waveguide. 

The lens shown in Fig. 1 initially focused the TMg mode at 92 mm as 
measured in the waveguide. After the 33 J/cm2 uv exposure, the focal 
length was reduced to 55 mm. The focal spot of the lens was scanned by 
coupling the transmitted beam out through a rutile prism and refocusing it 
onto an optical multichannel analyzer. The aperture was reduced to 3 mm 
in making the focal scans in order to ensure that all the light was 
captured by the relay lens. Traces of the OMA output are shown in Fig. 3. 

For a one-dimensional lens, the half-power diffraction-limited spot size 
is given by 
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( 1 ) 


S = 0.8859 O.Q f)/(next w); 

where A>q is the external wavelength of the light used, ngxt the mode 
index outside the lens region, taken as 2.289, f is the lens focal length 
as measured in the waveguide, and w is the input aperture. For the 
conditions stated, the expected minimTjm spot size after the uv exposure is 
4.5 pm. The measured value is 5.4 pm; so at this aperture the lens is 
about 1.2 times diffraction-limited. As the data in Fig. 3 show, before 
exposure to uv light, the focal spot was symmetric with most of the energy 
confined to the central peak, whose half-power width was 11.3 pm, 1.5 
times the diffraction-limited value, the first sidelobes being 6.7 dB 
down. After exposure to 33 J/cm2 of uv light the central spot width 
decreased to 5.4 pm, but some of the symmetry was lost and more energy was 
spread into the first sidelobes (now 5.0 dB down) and into the wings of 
the distribution. Further uv exposure shortened the focal length to 51 mm 
but led to deterioration in the spot quality. 

In summary, compact, fast, Luneburg-type lenses of good quality have 
been made by evaporation of AS2S3 onto LiNb03tTi waveguides through 
a rectangular mask, followed by exposure to ultraviolet light. There is 
room for improvement, though, in the focal-spot quality. 

We wish to thank Mark Parmenter, Richard Kenan, Carl Verber, and 
Philip Jenkins for their contributions to this research. Work at Battelle 
was supported by NASA-Langley Research Center and also benefited from 
support of parallel work by Harry Diamond Laboratories, U. S. Department 
of the Army. 
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Figure Captions 

Fig. 1. Twyman-Green interferogram of AS2S3 Luneburg lens on LiNb03 
substrate. The light source is the 488 nm line of an Ar laser, polarized 
predominantly vertically. The small interference features along the ridge 
probably result from diffraction due to the steep front and back edges of 
the lens . 


Fig, 2 , Focusing characteristics of a lens similar to the lens of Fig, 1, 
The transmitted beam was rendered visible by a white spray paint layer, 

(a) Normally incident beam, 9 mm wide, (b) Beam incident at 6.4 degrees 
from normal incidence, again 9 mm wide. The light source is the 633 nm 
line of a He-Ne laser, coupled into the waveguide with a rutile prism. The 
TE-polarized beam is propagating in the X direction in the crystal. 

Fig. 3, Photographs of CMA output, at input aperture 3 mm, showing 
intensity distribution in focal plane for lens of Fig. 1, Photos (a) and 

(b) show condensed and expanded horizontal scales, respectively, of the 
unexposed lens, while photos (c) and (d) show the condensed and expanded 
scales for the lens after exposure to 33 J/cm2 uv radiation. 
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1. INTRODUCTION 


The following report represents background research for, and technical 
embellishments of, a comprehensive computational model of a multispectral 
remote sensing system. This model is being developed as a tool for designing 
smart sensors which are able to process, edit and classify the data that they 
acquired. This report provides a forum for details and analyses that have 
not been suitable for reporting in other formats during the progress of the 
work. Consequently the reader will be assumed to be familiar with both the 
computational model and experimental results reported up to this date; (Huck, 
et al. , 1982 and Aherron, et al., 1981). In addition to the two above 
reports the AIAA (Breckenridge , 1979) and the SPIE (Barbe, 1979)’ have both 
published a collection of papers concerning smart sensors and the philosophy 
behind them. 

Accurately predicting the signal produced by a sensor observing a par- 
ticular target under a specific set of conditions is a very important goal 
for a model of this nature. Still more important is formulating a model 
that properly characterizes the stochastic properties of the signal that will 
likely be encountered. It is the physical nature of the interacting media 
that is considered to be variable and thus driving the derived quantities such 
as radiances. A quantity such as path radiance is not considered to have any 
inherent variability. Accounting for these variabilities introduced a special 
set of evaluation criteria in choosing a method of dealing with two of the 
major elements in the remote sensing system model: 1) atmospheric radiative 

transfer; 2) surface reflectance. 

In particular the implementation of the simulation required that the 
sensor irradiance be calculated for every simulated pixel with its associated 
random atmospheric conditions, or equivalently the equation of transfer 

t 

needed to be solved for every pixel simulated. This need of course discrim- 
inated heavily against lengthy iterative solutions to the transfer equation, 
accuracy not withstanding. An analysis of various atmospheric radiative 
transfer models will be presented. 

It is impossible txD over emphasize the importance of the spectral re- 
flectance of the targets in remote sensing work. Nevertheless, given the 



amount of research about and gathering of spectral reflectances taking 
place, particularly with respect to vegetation, it is amazing to find there 
is virtually no data report concerning the variability of spectral reflec- 
tances. Theoretical work concerning spectral reflectance has also been 
rather sparse. Several functional forms for reflectance variability will 
be examined and a collection of target reflectances will be cataloged. 

A goal established early in the development of smart sensor systems 
was the ability to detect cloud pixels. This ability would allow several 
options concerning data load reduction or choosing alternate imaging sites. 
The category of snow/clouds can usually be distinguished in either the visible 
or near- infrared. To distinguish between snow and clouds alone requires 
a spectral channel at a wavelength longer than 1.0 ym. 

The papers by Huck, et al. , (1982) and Aherron, et al., (1981) cover 
two sets of experiments performed over the period of this contract. The 
Huck, et al. paper concentrated on classifying pixels into fairly specific 
classes such as wheat and dark loam. The Aherron, et al. paper concentrated 
on the task of assigning pixels from any number of different substances into 
four broad categories which were: bareland, water, vegetation and snow/cloud. 

An enhancement implementing aerosol attenuation coefficients has produced 
results updated from those in Aherron, et al. A partial set of updated 
results for that work are presented. Results for a simple cloud detection 
algorithm will also be presented. 
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2. COMPUTATIONAL ATMOSPHERIC RADIATIVE TRANSFER 
2.1 Background 

The atmosphere is a significant element in the overall remote sensing 
model. A comprehensive treatment of this element Incurs the greatest com- 
putational costs of any model element. Numerous assumptions and simplifi- 
cations are usually made to facilitate ’’reasonable” limits to the complexity 
and cost of radiative transfer modeling. Different sets of these assumptions 
are manifest in the various atmospheric models currently available. Below, 
a brief discussion is presented of the problem and several computational 
models in use in remote sensing are highlighted and compared. From these 
discussions the justification for the chosen model should be apparent. 

Those not familiar with the equations and terminology of atmospheric radia- 
tive transfer should consult Liou (1980) or Wolfe and Zissis (1978) for a 
background development. The terminology and notation used in this and 
related papers is almost identical to that used in Chapter 4 of Wolfe and 
Zissis (1978) which makes it particularly good as an Introduction. 

The fundamental equation of remote sensing after Slater (1980) 

L = i E T U p T + L p T + L (2-1) 

71 o o o a p 

where E^ = is the solar spectral irradiance at the top of the atmo- 

sphere, T^ = T^(X,t,P^) is the atmospheric transmittance along the incident 
path from the sun to the surface (solar zenith angle “ ^^9 ®o^ ^ 

E L^(E^,X,T,y^,p) is the diffuse sky spectral radiance which results from 
all radiation scattered downward onto the surface (i.e., integrated at the 
target over elevation and azimuth); p = p(X) is the spectral reflectance 
of the surface (sometimes called ’’signature”); T = T(X,T,y) is the atmosphere 
transmittance along the exitant path from the surface to the sensor (zenith 
angle = 0, y = cos0) ; and L = Lp(E^,X,T,y^,y,(t)) is the path spectral rad- 
iance which results from all radiation scattered upward along the path from 
the surface to the sensor. The other parameters are wavelength, X, optical 
thickness of the atmosphere t = x(X), and azimuth angle (|) between the planes 
of incidence and exitance. The component of the total radiance L which arises 
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from radiation reflected from the target is referred to as the beam spectral 
radiance L. , that is L, = L - L , and L, = L, (E ,X,T,y ,y,p). 

D D ' p D D O O 

Now the equation of radiative transfer. for a plane parallel atmosphere 
(Malila, et al. , 1971) is given by 

,T w (t) „ 

y^ = L(T,y,(|)) rVJ p(T,y,<j),y',({>’)L(T,y*,<j)')dy’d<i)’ - (l-w^(T))B(T) 

( 2 - 2 ) 

where is the single scattering albedo which is the ratio of the sum 

of the scattering coefficients to the sum of the scattering and absorption 
coefficients. B(t) is the Plank radiation function for thermal self-emission. 
The single-scattering phase function p(*) can be described as a probability 
density function for the particular direction a photon will scatter relative 
to the original direction of travel. The function usually has a highly 
irregular (radially asymmetric) shape. The straight numerical solution of 
this equation is very time consuming. The first assumption made for even 
the complex form shown is that the atmosphere is plane-parallel, that is it 
is an infinite slab bounded below by the ground. For near-nadir looking 
sensor systems this is not an unreasonable assumption. Several uses of 
atmospheric radiative transfer models in the remote sensing literature and 
the methods of a solution for equation will now be examined. First, atten- 
tion will be focused on the solution of the equation of radiative transfer 
as implemented in the various methods. Assumptions and approximations con- 
cerning the atmosphere will be discussed later. 

2.2 Radiative Transfer Models 

The models discussed will be limited to those that have received 
attention in the remote sensing literature. Probably the most widely known 
model is the one developed by Turner (Malila, etal., 1971 and Turner, 1974) used 
in early systematic remote sensing studies. Another model receiving more 
recent attention was developed by Dave (1978) and has been applied to both 
remote sensing problems, Dave (1979), and solar insolation. O’Neill, et al. 
(1977, 1978) implemented Liou’s (1973) model and compared it to Turner’s 
model and to actual LANDSAT data. Kiang (1980) implemented Hansen’s (1969) 
model for tests concerning atmospheric effects on Thematic Mapper data. The 
other less general models were encountered, due to Otterman, et al. (1980), and 
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Lampley and Blattner (1978), upon which discussion will be delayed. Attention 
will be concentrated on those models that account for all the terms in Eq. 
2-1. In dealing with the atmosphere in strictly a transmissive mode, that 
is accounting for only the first term in Eq. 2-1, LOVTTRAN (Selby, et al. , 
1978) has become the accepted standard model. 

The fundamental character of, . and difference between most approximate 
methods of solution for the full equation is in the form and order of 
approximation of the phase function. This approximation is to allow the 
elimination of the integration in Eq. 2-2. The phase function determines the 
directional characteristics of all scattered radiation. Wolfe and Zissis 
(1978) show plots of some typical phase functions that demonstrate their 
pronounced asymmetry. Turner used two Dirac delta-functions one forward 
and one backward to approximate the phase function. Liou (1973) expanded 
the phase function into a series of N Legendre polynomials in his Discrete 
Ordinate Method (DOM). References to this model frequently are of the form 
**2N-stream approximation** where N is an integer. Dave (Dave and Canosa, 

1974) used the method of Spherical Harmonics where the phase function and 
the intensity were expanded into a series of L + 1 Legendre polynomials. 

Krook (1955) shows that the DOM and Spherical Harmonic methods are equivalent. 
The accuracy of these equivalent methods can be made arbitrarily high by 
increasing the number of terms in the expansion (O’Neill, et al. , 1977). 

The typical scattering phase function is a highly **peaked** function and 
thus rich in harmonics. In essence the Legendre based methods approximated 
the phase function with some finite set of polynomials (i.e., harmonics) 
while the Turner method approximated the function with an infinite set of 
harmonics in the form of the dirac-delta function. All of the above methods 
can be defined for a single homogeneous layer of any reasonable optical 
thickness. Hansen (1969) approximated the phase function at 20 discrete 
directions. The other f undemental difference is the intermediate solutions 
reached are valid only for an optically thin layer where multiple scattering 
can be ignored. The properties of an optically thick layer are arrived by 
aggregating the thin layers. Hansen’s model accounts for polarization which 
is unique for the models discussed. 
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2. 3 Model Comparisons 

Ronnholm, et al* (1980) implemented a form closely related to each of 
the models discussed. The implementation of the radiative transfer irwdels 
by Ronnholm, et al. assumed several things about the atmosphere. First an 
analytical form for the "real” phase function was used (Henyey-Greenstein) . 
Second, the atmosphere was assumed to be vertically and horizontally homo- 
geneous. All of the models except Turner’s have been implemented with 
vertical inhoraogeneity capabilities but this restriction was necessary for 
comparison. Each model was run for the same set of atmospheric conditions. 
(The Ronnholm, et al. literature citations will not be duplicated). Ronnholm, 
et al. used Twomey’s doubling method model as a benchmark or "truth" so to 
speak, considering it to be more accurate than the other methods compared. 
Hansen’s (1969) and Twomey’s models are comparable. A delta-Eddington 
approximation by Joseph was included which is similar to the Turner model. 

The Discrete Ordinate Method for 4 streams was implemented using the ana- 
lytical solution (Liou, 1974). Since the form of the differential equations 
solved for the DOM and Spherical Harmonics methods are similiar, the DOM’s 
computational times will be considered representative of the Spherical 
Harmonics Method (when L = N) . 

One of the major findings of Ronnholm, et al. was the relative compu- 
tational burden for the three classes of methods. The results are sum- 
marized in Table 2.1. The results in Table 2.1 can only be meaningful 
within the context of the relative accuracy of the various methods. Ronnholm, 
et al. compare the three methods for various conditions and concludes "... the 
factor of 20 in computation time saved by the delta-Eddington was judged more 
valuable than the greater accuracy of the four-stream method." They add 
"... if optical depths, single scattering albedos or asymmetry factors are 
either uncertain or known to be flucuating;. with relative standard deviations 
of 10% or greater, then little real benefit is added by the use of compu- 
tationally precise, but costly, many-stream radiation-transfer algorithms." 

An important point that should be made here, is that Ronnholm, et al. com-- 
pare values of reflection, diffuse transmission, and absorption for the 
atmosphere and thus net fluxes are being compared. For aerosol laden atmo- 
spheres, and their highly anisotropic scattering properties, the jump from 
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TABLE 2.1 


Relative computation burden of radiative 
transfer models used by remote sensing 
investigators analogy to methods implemented 
by Ronnholm, et al. 








comparing fluxes to comparing intensities with directional dependences is 
not straight forward. For example, the radiation field incident on a scat- 
tering volume has only two components in a two stream approximation and for 
highly anistropic scattering this is an inadequate approximation for calcu- 
lating the scattering into the highly specific direction of the sensor. In 
general, the more anisotropic the scattering phase function, the more accu- 
rately it needs to be approximated in order to determine the radiation field 
within the atmosphere. For the Turner model, the radiation field within the 
atmosphere is approximated utilizing the double-delta phase function. The 
actual path radiance contribution of a scattering volume is calculated using 
that simplified radiation field and interpolated values of the scattering 
phase function. The hemispherical integration for flux quantities implicit 
in Ronnholm, et al. tends to mask the inaccuracies. This provides some 
explanation for the disagreement between Ronnholm, et al., who found little 
difference between methods for flux quantities, and O^Neill, et al. , who 
found significant differences between the predicted path radiances for the 
DOM and Turner’s model. Even so the remarks by Ronnholm, et al. concerning 
the uncertainty about atmospheric conditions should not go unheeded, par- 
ticularly for the case of errors in approximations of the scattering phase 
function. 

It is useful to elaborate on the efforts of O’Neill, et al. concerning 
calculation of path radiance (O’Neill, et al. 1977, 1978). In O’Neill, 
et al. (1977) Turner’s method and Liou’s DOM were implemented assuming 
vertical homogeneity. The results were fairly comparable for low aerosol 
optical depths but were divergent as optical depth increased. They attrib- 
uted this to the crude angular approximations used by Turner and their 
increasing importance at greater optical depths. In O’Neill, et al. (1978) 
path radiances computed by the DOM were found to be in good agreement with 
path radiances determined from LANDSAT data using clear lake reflectors. 

The DOM performed considerably better than the Turner method for the same 
set. of data. Horvath, et al. (1972) found that the Turner method had 
systematic errors in predicting atmospheric effects on aircraft multispectral 
scanner data. The input data used in Horvath, et al. was not of as high a 
quality as for the O’Neill, et al. work. 
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The assumption of vertical homogeneity within the atmosphere does not 
seem plausible considering the distribution of aerosols. Neither Ronnholm, 
et al. (1980) nor Liou (1975) could demonstrate much difference in modeling 
results, for flux related quantities, between the lesser and more "accurate” 
methods of calculation. No test of the assumption's effect, in terms of 
model results, on radiance values was found. O'Neill, et al. achieved good 
agreement between LANDSAT determined and model predicted path radiances 
under the vertical homogeneity assumption. 

By far the more important atmospheric constituent in terms of its 
effect on remote sensing is the aerosols. In certain wavelength regions 
though, absorption by other constituents becomes important. The major 
ones being water vapor, ozone,, and oxygen in the visible and near infrared 
regions. Accurately accounting for absorption by the constituents is not 
a trivial task (Selby, et al., 1978; Dave, 1978). Turner treats the ozone 
absorption, occuring in the upper atmosphere, as a phenomena separate from 
the scattering, which takes place in the lower atmosphere. Except for Turner 
those models that do deal with molecular absorption utilize the Air Force 
Geophysics Lab's model LOWTRAN (Selby, et al. , 1978) or data from that model. 

2.4 Special Models 

As was mentioned before, two other less general models of transfer 
were found in the literature. The Otterman, et al. (1980) work was mostly 
illustrative in nature, dealing only with low optical thickness atmospheres, 
and is rather unique in its solution. Both Otterman, et al. and Lampley 
and Blattner (1978) deal only with single scattering. This assumption 
greatly simplifies the characterization of the radiation field within the 
atmosphere. For the case of optically thin atmospheres this is not a bad 
assumption but likely to be encountered only in more arid regions on a 
regular basis. This assumption, for thin layers, is the basis of the dou- 
bling method discussed earlier. Eliminating the effects of multiple scat- 
tering is. .judged, to-be. too large a source of error when dealing with atmo- 
spheres optically thick enough to cause significant image degredation. 

2.5 Model Selection 

One of the evaluation criteria set for evaluation of various atmospheric 
models was their suitability for use as a computer subroutine for a remote 
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sensing system simulation. One other approach suggested by O'Neill, et al. 
(1978) and inferred by Dave (1978) is the use. of table-lookup method. For 
certain problems this could be a very useful technique. But for a general 
purpose source of radiance data for a stochastic simulation, with many of 
the input parameters being random variables, it would likely prove excessive 
in its storage requirements. To define the size of radiance data base needed 
to satisfy the modeling requirements would require knowledge as to how fine 
a grid, or subdivision, in each parameter, would be necessary to retain the 
desired accuracy. This knowledge is not available a priori and would only 
come from careful and laborious sensitivity analysis for the input para- 
meters, It may be possible to gain increased accuracy, over say the Turner 
model, by implementing a table lookup scheme for a more accurate model 
but it would be better to pay the greater price for the subroutine method 
and get the full benefit of increased accuracy. 

From what has been discussed one can see that the Turner model of 
radiative transfer introduces significant errors to the approximation of 
sensor incident radiation compared to other more comprehensive models (i,e. 
DOM). There is though, an estimated twenty fold reduction in computing 
time realized, I^ether or not the error introduced by the Turner model is 
significant in the light of the uncertainties in the atmospheric parameters 
has yet to be determined. Ronnholm, et al. showed that, for flux calcul- 
ations, these errors were not significant. For an initial Implementation 
of a model with multiple scattering capabilities within a comprehensive 
remote sensing system model, the Turner model is the logical choice- Once 
greater knowledge of the interrelationships and sensitivities of the various 
system entities is established, it would probably be best to upgrade to 
the Dave model. 

2. 6 Model Implementation 

The Turner model was implemented as a subroutine within the computer 
simulation code. The model is available for purchase from the Environmental 
Research Institute of. Michigan (ERIM) located in Ann Arbor, Michigan. Scat- 
tering phase function data files which are distributed with the model were 
used. Input to the Turner model, as implemented, consists of nine parameters 
and one function (table form). Those parameters considered to be deterministic 
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(fixed for a particular run) were - Solar Irradiance, 0^ = Solar Zenith 
Angle, (}) = Azimuth angle between incident and exitant vectors, 0 = Sensor 
Zenith Angle, 0)^ - Single Scattering Albedo. The scattering phase function 
was also deterministic. All simulation to date have been performed with 
phase functions from a continental type aerosol model with a complex index 
of refraction of m = 1.5 - .Oli. The single-scattering albedo was chosen 
to be .9 and is constant. Simulations have been confined to nadir looking 
sensors with solar zenith angles in the range of 30° - 50°. 

There are four stochastic parameters consisting of = background re- 
flectance, = target reflectance, = Rayleigh optical thickness and 
= total optical thickness. The background used was a bare land target. 
Both reflectances are generated according to the functional form discussed 
earlier. The optical thicknesses are generated by the method discussed in 
Buck, et al. (1982). It is necessary to separate the Rayleigh component 
because the phase function used is; a weighted mean of the Rayleigh and 
aerosol phase functions. 

The single-scattering albedo used is constant over wavelength which 
requires a word of warning. ' The single-scattering albedo defines the rela- 
tive mix of scattering and absorption in the contribution to attenuation. 
Therefore in moving from spectral regions characterized by scattering to 
those with significant absorption the single-scattering albedo changes. 

A method for varying the single-scattering albedo has not been implemented^ 
The single-scattering albedo chosen (.9) is reasonable for regions of the 
spectrum without significant absorption. It is therefore necessary to 
restrict radiative transfer modeling to non-absorbing regions at the present 
time. 
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3. SPECTRAL REFLECTANCE MODELS 

Ultimately it is the spectral reflectance characteristics of targets, 
perhaps in concert with spatial distribution characteristics, that provides 
the information to users of remotely sensed data. It is therefore important 
that the target reflectance properties introduced to a simulation be as 
realistic as possible. This is particularly true for optimization studies 
where optimizing a system for a set of "artificial" targets is of more 
academic interest than practical use. 

There were two major tasks to be addressed in the area of reflectances. 
One concerned characterizing the deterministic nature of spectral reflec- 
tances which consisted of assembling representative spectral reflectance 
curves for a number of targets of interest. The other task dealt with 
characterizing the stochastic nature of spectral reflectances by selecting 
and parameterizing functional forms for describing reflectance variability. 

3. 1 Reflectance Functions 

The quantity utilized in the simulation is. the spectral diffuse re- 
flectance defined by Slater (1980). 

■ E (X) 

where the target is assumed to have Lambertian characteristics. The quan- 
tity actually measured in the field, depending on the particular measurement 
program, is more likely the hemispherical-conical reflectance factor which 
is then normalized to some standard reflector (Smith and Ranson, 1979). 

Present work has been limited to the Lambertian target assumption though 
natural targets have been demonstrated to deviate significantly from this 
assumption (Smith and Ranson, 1979). In general Slater (1980) offers a 
very good introduction to surface reflectance and its importance in remote 
sensing and presents a brief introduction to several efforts aimed at 
modeling plant canopy reflectance. Smith and Ranson (1979) offers a fairly 
comprehensive review of the literature concerning data and/or models of the 
directional reflectance characteristics of natural surfaces. 

In order to produce spectral reflectances with stochastic characteristics 


X 
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(3-1) 


it was necessary to choose an analytical form for reflectance. The func- 
tional form chosen for the reflectance is given by 


-B (A)X 

p(X) = p (X) e ° 
o 

This is the same functional form (Bouguer's Law) as used for atmospheric 
transmission. The wavelength dependence of reflectance will be omitted in 
subsequent notation. This form was used because of its analytical tract- 
ability and compatibility with the concurrent atmospheric modeling efforts 
and because the family of curves generated were qualitatively similar to 
sets of field data (Park, et al., 1980). The parameters of the distribution 
are P^(^) and B^(X). This is identical to a form for spectral reflectance 
proposed by Tucker and Maxwell (1976) for vegetation canopies in regions 
of the spectrum characterized by canopy absorption (i.e. , the chlorophyll 
absorption regions). Within the same spectral region and dependent on the 
value of X, Tucker (1977) and Tucker and Maxwell (1976) found that 


P - 


(3^2) 


was also a useful functional relationship. 

In Tucker’s work X represented any one of several plant canopy char- 
acteristics including several biomass measures, leaf water content and 
chlorophyll content. For the near infrared region from .74 m to 100 m which 
is characterized by little absorption and relatively high reflectance due 
to scattering. Tucker (1977) found, again dependent on the value of X, the 
following two equations to be good approximations of reflectance 


p = Po + 


(3-3) 


p = S 


(1 


Po + 

- . e 


6 X, , 

O) 


(3-4) 


Where S represents the asymptotic reflectance of the vegetative canopy. 

This is the reflectance that an Infinitely thick canopy would display. Park 
and Deering (1981) developed a set of differential equations (modifications 
of the Kubelka-Munk model) for describing the interaction of light with 
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plant canopies that is remeniscent of the Eddington Approximation for atmo- 
spheric radiative transfer. If the background reflectance is taken to be 
zero then the diffuse reflectance formula of Park (retaining previous 
notation) 


P = 



-3 X 
o o 



X 


o o 


(3-5) 


Of the above functional forms, only Eqs. 3-2, 3-3, and 3-4 display 
asymptotic nature. Equation 3-2 asymptotes for large values of X and 
Eqs. 3-3 and 3-4 asymptote for small and large positive values of X. 

Ideally if X is some vegetative measure such as biomass per unit area then 
the reflectance functions should asymptote to the background reflectance 
for small X and approach the vegetative reflectance asymptote for large X. 
These three equations were formulated to be driven by inherently positive 
plant characteristics such as canopy weight. Equation 3-1 is driven by 
an artificial variable taking positive and negative values. As was pre- 
viously stated, the resulting curves were qualitatively acceptable. 

Target reflectance phenomenon represent a special case of radiative 
transfer, a subject which has received a great deal of attention over the 
years. Commenting on the dearth of surface reflectance modeling efforts 
Smith and Ranson (1979) hypothesize "probably, this is a recognition of the 
difficulty of specifying the appropriate phase function in both a sufficient 
and tractable manner and further, performing the necessary measurements to 
determining the phase functions". They were commenting on the difficulty 
of specifying and measuring the phase function of such things as twigs and 
leaves. Smith and Ranson (1979) discuss the formulation development of two 
major vegetative canopy models, Smith and Oliver (1972) and Suits (1972), 
taking care to link the "... analytical and physical reasoning of canopy 
radiation interactions to the broader mainstream of radiative transfer 
theory". A better understanding and overview can be gained by reading 
their unified discussion prior to working with the original papers. The 
radiative transfer methods used in these models can be classed as those that 
use simplistic approximations to the phase function such as the Eddington 
approximation. 
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The following development outlines the parameterization of Eq. 3-1- 

We assume X is a (standard) normal random variable with mean u =0 
o p 

and variance =.l- Repeating the formula for reflectance- The assumption 
of normality for X will be discussed more completely in a later section. 

When X is assumed to be a normal variate the reflectance, given by Eq. 3-1, 
has a log-normal distribution- 

0 ' Po' (3-6) 


then 





1 ) 


(3-7) 

(3-8) 


If we have estimates for y and a denoted y and a , respectively, 

P P P P ^ ^ 

then by substituting the estimates into the above formulas yields estimates 

of the parameters. 




(3-10) 


Finding data to define the variability of spectral reflectances (i.e. , 
Op) has not been very fruitful. As of yet, no data has been found, on the 
covariance (between wavelengths) of spectral reflectance for various targets - 
Values for the variance (i.e., the diagonal terms of the covariance matrix) 
can be estimated from data, for wheat, reported by Collins (1978) through 
the following reasoning as shown below. 

A simple expression for the remote sensing equation for vertical sun 
and sensor can be defined as 
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L = ^ E p + Lp 

(With notation as before in atmospheric section). Over small geographical 
areas it can be assumed that all the terms are constant except p. Allowing 
the calculation of 

a = - E a 

L TT O p 
L = - E p + Lp 

Collins reports data for medium altitude aircraft in a fairly dry area 
(Imperial Valley) thus Lp is ignored and thus dividing Eq- 3-12 by Eq. 

3-3 gives 


Ol a 

= Coefficient of variation (CV) 

L p 

Collins reports CV as percent standard deviation of measure radiance which 

Op 

determines necessary for formulas 9 and 10. Since Collins was working 
with radiances an estimate of p necessary for Eq. 3-10 is not available. 

Rao, et al. (1978) similarly reported reflectance CV’s for grain crops 
and soil. The data was for variation over several months and encompassed 
atmospheric corrections of unknown nature making it unsuitable for present 
efforts. A CV of .1 was chosen for all agricultrual crops after reviewing 
Collins data. 

As a note of explanation, Smith and Oliver (1972) developed a stochas- 
tic vegetation canopy model but more appropriately should be termed a Monte 
Carlo model in that the photon interaction with various canopy constituents 
(canopy orientation and distribution) was treated as a random process but 
then the canopy composition was constant and thus is fundamentally different 
than the model being discussed. 

So far discussions of reflectance have centered on vegetation. The 
same function was also used for the other targets. For soil Condit (1970) 
was a valuable source of spectral reflectances. The nature of the curves 
presented allowed an estimate of reflectance variability to be made. Each 
soil was characterized by two spectral reflectance curves: one for dry 
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soil and one for wet soil. The mean and variance for that soil was deter- 
mined by the following . 




"dry 


wet 


Pdrv ~ Pwet 


In general, variances were chosen for other targets in order to achieve 

"Reasonable families" of curves in a stochastic simulation. One significant 

modification of the parameters was made for vegetation. The reflectance 

variability was made proportional- to reflectance. That is CV was constant 

(.1) and multiplied by p to yield for estimating 3^ and p^. The effect 

of this modification can be seen in the difference in variability behaviour 

between vegetation and sand in Fig. 4 of Huck, et al. (1982). ■ As one can 

see this introduces variations in O over wavelength. It should be noted 

P 

that the ± sigma plots of Fig. 4 (Huck, et al., 1982) are for ±a and not 

±a such that they enclose -67% of the simulated reflectance values even 
P 

though p is not normally distributed. This should become clear in the 
section on statistical characteristics. 

The reflectance covariance matrix for targets in general is defined 
as 


1=01 
P P n 


for 


i = 1, 2, 3, . . . , n 


and where .1^ = n x n identity .matrix -and^ 0^. is constant, -but defining 


-2 
P = 


-2 


0 


'1 -2 


'2 -2 


_ 9 

P^ 


n 
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a 

and CV = ^ 

P 

the covariance matrix for vegetation is 

2-2 
E = CV I p 
p n 

where CV is taken to be constant. 

It is particularly important that the off-diagonal terms of the reflec- 
tance covariance matrix be approximated in future work and especially for 
investigation concerning the optimum choice of spectral channels. It has 
frequently been observed, or bemoaned if you will, that several of the 
LANDSAT channels are redundant. This is due in part to the high correlation 
of reflectance in the different channels. For example both LANDSAT channel 

6 (.7 - .8) and channel 7 (.8 - .11) fall in the spectral region for vege- 
tation characterized by the same reflectance phenomenon. Thus changes in 
the nature of the plant. that affect reflectance in channel 6 affect channel 

7 in almost the same way. This spectral reflectance correlation for channel 
6 and channel 7 is true for most targets and therefore one of the channels 
is fairly redundant. 

Wiersma and Landgrebe (1979) make an important contribution to remote 
sensing where they “develop an analytical procedure for the design of the 
spectral channels for multispectral remote sensor systems". Though never 
explicitly stated, it appears their analysis is based on spectral radiance 
data. The important point to this discussion is that they chose an optimal 
(mean-square) set of spectral channels where the analysis was based on the 
spectral radiance covariance matrix and is aimed at reducing channel re- 
dundancy. By reducing correlation between the channels they increase the 
information available from a fixed number of channels. This is shown by 
the techniques of Information Theory. Their radiance data does not allow 
the estimation of the spectral reflectance covariance necessary for the 
present modeling efforts. 

3 . ] Target Reflectance Curves 

As was indicated, the target reflectances represent an important el- 
ement in the remote sensing simulation. A collection of target reflectances 
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was gathered and subsequently used in simulations reported to date. Data 
collected was limited to in situ reflectances for single or collections 
of whole plants. Only data available in the open literature was utilized. 
There exists a large body of reflectance data produced for NASA, in raw 
form, that has yet to be utilized. 

Examination of remote sensing literature from the early 60* s to present 
shows the tendency to report reflectances to longer wavelengths in the in- 
frared. Even with this tendency, very little reflectance data is published 
for wavelengths longer than 1.1 ym which happens to correspond to the limit 
of LANDSAT spectral coverage. Some data utilized in the cloud detection 
tests, which required data out to 1.60 ym, was of rather crude spectral 
resolution especially in the .50 ym to .80 ym range that is so important 
for the BAM categorization discussed in Aherron, et al. (1981). The sources 
of reflectance data are listed in Table 3.1 along with the figure number (s) 
for the plot(s) of that data. 
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TABLE 3.1 
REFLECTANCE DATA 


FIGURE NUMBERS 


SOURCE OF DATA 


FIG. 

3.1 - 

3.4 

Lansing, 1970 

FIG. 

3.5 - 

3.6 

Condit, 1970 

FIG. 

3.7 


Vlcek, 1974 

FIG. 

3.8 


Suits and Safir 


FIG. 3.9 (Silt Loam) Bowers and Hanks, 1965 


FIG. 3.9 
FIG. 3.10 
FIG. 3.11 
FIG. 3.12 
FIG. 3.13 


Orr, et al. , 1963 
O'Brien and Munis, 1975 
Hansen, 1969 
Novosel' tsev, 1965 
Leeman, et al. , 1971 






REFLECTANCE 


WHEAT HIGH 
WHEAT LOW 


OATS I 

OATS It 



WAVELENGTH, 


FIGURE 3.1 SPECTRAL REFLECTANCES 


1.0 

.8 


BARLEY H 

BARLEY M 

TOMATOES 

PINE 

POPLAR 


.6 



WAVELENGTH, firu 


FIGURE 3 . 2 SPECTRAL REFLECTANCES 
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REFLECTANCE REFLECTANCE 



.3 ,5 .7 .9 1.1 1.3 1;5 1.7 1.9 2.1 2.3 


WAVELENGTH, 

figure 3.5 SPECTRAL REFLECTANCES 



WAVELENGTH, /im 

FIGURE 3.6 SPECTRAL REFLECTANCES 
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RED PINE 
WHITE PINE 


/ 
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.9 1.1 1.3 1.5 1.7 1.9 

WAVELENGTH, [im 



2.1 2.3 


FIGURE 3.7 SPECTRAL REFLECTANCES 
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FIGURE 3.9 SPECTRAL REFLECTANCES 
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FIGURE 3.10 SPECTRAL REFLECTANCES 
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REFLECTANCE 



WAVELENGTH, jim 


FIGURE 3.13 SPECTRAL REFLECTANCES 
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4. CLOUD DETECTION 


4,1 Background 

For an Earth resource observation system clouds present a major 
limitation to the collection of data. Worldwide cloud cover is approximately 
45% (Barnes et al., 1968), which represents a very significant potential 
for obscuration of ground targets. The ability to indentify clouds reliably 
would allow several options for smart sensor systems including the ability 
to avoid imaging of cloudy areas and to eliminate cloud pixels from the data 
to be transmitted to ground stations. An even more sophisticated task for a 
sensor system with reliable cloud detection capability would be the selection 
of alternate imaging areas in the event that the primary target is cloud 
covered . 

In order to identify clouds they must be adequately characterized spect- 
rally. As one can imagine, determining the spectral reflectance of a cloud 
is no easy task. The reflectance data available for clouds has generally 
been the product of radiative transfer modeling (Hansen 1969, Novoseltsev, 
1965). In the visible and near-infrared region of the spectrum both snow 
and clouds form a single category of what might be called bright white 
targets. Within this spectral region the two classes have such closely 
similar spectral characteristics as to make the ' identification of clouds 
as a class unreliable. Only in the longer infrared wavelengths do the two 
classes show significant differences. 

Valovcin (1978) made a important comprehensive study of the snow/cloud 
discrimination problem specifically aimed at potential on-board methods. 
Valovcin utilized, for his analysis, spectral radiance data from cummulus 
and cirrus clouds and snow obtained from a high flying aircraft. He analyzed 
several different discrimination tasks including identifying different cloud 
types. For the simple cloud versus snow discrimination he found that a 
radiance threshold for a .11 ym bandwidth spectral channel centered at 1.56 ym 
achieved an 85% + accuracy. Thoygh no. test results were reported for such, he 
suggests that moving the channel center to 1.525 ym might improve results. 

4.2 Test 

Based on results and data reported by Valovcin (1978) a new channel was 
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added to the basic set of two channels used in the target categorization 
tests reported in Aherron, et al. (1981) • Those channels consisted of one 
centered at .65 ym and one centered at .85 ym with a bandwidth of .02 ym. 

The new channel for cloud detection was centered at 1.55.yniwith a . 1 ym 
bandwidth. 

As was indicated above the reflectance data utilized in the modeling 
was the result of theoretical investigations. Hansen (1969) used the doubling 
method solution to the radiative transfer equations, discussed in the section 
concerning atmospheric radiative transfer, to determine the diffuse spectral 
reflectance of ice clouds of various optical thickness. Novolsel* tsev cal- 
culated the spectral reflectance of a cloud composed of water droplets. 

Zander (1966) measured the reflectance of laboratory generated ice clouds and 
his data compares favorably with that report by Hansen (1969). Stochastically 
clouds are treated the same as non-vegetated surfaces. 

A cloud detection algorithm was formulated based on the 3 above channels. 
The algorithm is basically a two level decision tree. First, the snow/cloud 
category is separated from the other categories based on a threshold radiance 
in the . 65 ym channel. This is the method previously used in simple categori- 
zation (Aherron, et al. 1981). The second level decision of snow versus cloud 
is based on a threshold in the 1.55 channel. Figure 4.1 shows the covariance 
elipse plots for a collection of targets similar to those used in Aherron, et al. (1981) 
simulated for 23 km visual range. The first level categorization boundaries 
are indicated by dotted lines. Included in the snow/cloud category are a snow 
target, an optically thick water droplet cloud, and four ice clouds of various 
optical thicknesses. The ice clouds are indicated on the ellipse plot accord- 
ing to optical thickness. 

The ill-placed clouds are a manifestation of two problems in dealing with 
clouds. First is the problem of thin clouds. It is an arbitrary decision as 
to what constitutes a cloud and what is merely haze in terms of the magnitude 
of their effect on the received signal. Secondly, the method by which clouds 
are simulated assumes a perfectly absorbing background (immediately behind the 
target as opposed to that outside the field of view). Thus clouds of low 
optical thickness, with a corresponding significant transmissive contribution 
to radiance will be simulated as having too low a radiance , by a magnitude 
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CHRNNLL 1 ( .65 ym) 



FIGURE 4.1 COVARIANCE ELLIPSE PLOT FOR 23 km VISUAL RANGE AND 30% 
SOLAR ZENITH ANGLE WITH ICE CLOUDS ANNOTATED WITH 
OPTICAL THICKNESS. 
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dependent on the intended background. The cloud set used to test the algorithm 
was limited to those with optical thickness of 16 or greater in order to make 
the results of snow versus cloud more meaningful. This limits the test of 
the ability to distinguish snow versus clouds to those pixels that can be 
reliably categorized as snow/clouds. 

Figure 4.2 shows the ,65 pm versus 1 .55 yra ellipse plots for a similar but 
smaller set of data as used for Figure 4.1. The cloud and snow classes have 
been indicated. A threshold value in the 1 .55 yra channel can be seen to 
separate snow and clouds fairly well. This algorithm was implemented and 
tested for the three atmospheric conditions (as qualified by visual range), 
and used two solar zenith angles for testing in Aherron, et al. (1981). 

Overall the snow versus cloud discrimination accuracy was very good. 

The worst case was for a visual range of 5 km but the difference between 
5 km, 10 km, and 23 km results was insignificant. Table 4.1 shows the snow 
and cloud confusion matrix for the 5 km case. 

The above test has served to show that given perfect snow/cloud versus 
other categorization the snow versus cloud discrimination is simple and ac- 
curate. In practical applications the snow/cloud versus other categorization 
can be rather unreliable. It will be the major source of error in dealing 
with clouds. 
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SNOW 
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FIGURE 4.2 COVARIANCE ELLIPSE PLOT FOR 23 km VISUAL RANGE AND 30° 
SOLAR ZENITH ANGLE WITH SNOW AND CLOUDS INDICATED 
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TABLE A.l 


TRUTH 


Snow versus Cloud confusion matrix for 
5 km visual range and 30“ solar zenith 
angle. 

DECISION 



SNOW 

CLOUD 

SNOW ! 

198 

2 

CLOUD 

0 

346 




5. MODEL ENHANCEMENT TESTS 

Subsequent to the publication of Aherron, et al. (1981) an enhancement 
was made to the aerosol and Rayleigh attenuation coefficients (ai(A) ) which 
altered the spectral radiance signatures produced by the simulation. Figure 
4.1 represents signatures produced utilizing the new coefficients. Two of 
the three experiments in Aherron, et al. were repeated with enhanced coef- 
ficients and are summarized in Table 5.1. 

Table 5.1(A), which corresponds to Figure 8(B). of Aherron, et al., shows 
results for the 7 class set. There are two major conclusions that can be 
dravm from this figure. One is that. the MLH (Maximum Likelihood) aggregation 
method is fairly sensitive to changes in solar zenith angle. Second, recalling 
that the BAM (Boundary Approximation Method) is substantially cheaper computat- 
ionally, the BAM is the best method of the three for the task of categorization 
for this set of classes. 

Figure 5.1(B), which corresponds to Figure 8(C) of Aherron, et al. , 
shows results for a larger set of 20 classes. Again the MLH method proves 
susceptable to changes in solar zenith angle. The MSD (mean square distance) 
aggregation method gave the highest accuracy but with the increased cost over 
the BAM. These new results show even more dramatically the importance of the 
mix of classes ( i.e taisk assignment) in determining sensor system performance. 
Ultimately the mix of classes on which sensor systems are tested should be 
based on realistic probability of occurrence data. This final point will be 
elaborated on in Section 6. 
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. TABLE 5.1 (A) 

CATEGORIZATION ACCURACIES FOR 3 VISUAL RANGES 
AND TWO SOLAR ZENITH ANGLES FOR 17 CLASSES. 


SOLAR ZENITH ANGLE 


VISUAL 

RANGE 

30“ 

— 

O 

o 


MLH 

MSD 

BAM 

MLH 

MSD 

BAM 

23 km 

1. 

.96 

.99 

.91 

.96 

.98 

10 

1. 

.96 

.99 

.91 

.95 

.98 

5 

1. 

.96 

.99 

.91 

.96 

.98 


TABLE 5.1 (B) 

CATEGORIZATION ACCURACIES FOR 3 VISUAL RANGES 
AND TWO SOLAR ZENITH ANGLES FOR 20 CLASSES. 


SOLAR ZENITH ANGLE 


VISUAL 

RANGE 

30“ 

o 

o 


MLH 

MSD 

BAM 

MLH 

MSD 

BAM 

23 km 

.96 

.93 

.86 

.85 

.90 

.86 

10 

.95 

.91 

.86 

.85 

.90 

00 

• 

5 

.96 

.92 

.87 

.85 

.91 

.86 


ALL TRAINING SETS GENERATED FOR 23 km VISUAL RANGE 
AND 30“ SOLAR ZENITH ANGLE. 





6. STOCHASTIC MODELING 


6,1 Introduction 

This section deals with the stochastic nature of the remote sensing 
system simulation. The analytical details of the stochstic properties of 
the model have been discussed in other sections or in other papers (Huck, 
et al. , 1982). Included in this section are observations and discussions 
concerning the limitations encountered and the implications of various 
assumptions made in the treatment of stochastic model properties. 

The general philosophy in addressing the statistical properties of 
the remote sensing system has been to statistically characterize the fund- 
amental values or processes within each model element. For instance, during 
the discussion of atmospheric radiative transfer models, it was stated that 
it was the amount and characteristics of the atmospheric constituents that 
was “driving” variation in radiance quantities. The limit of resolution to 
which this fundamental characterization is carried is defined by practical 
realization, as opposed to availability considerations. For plants it is 
realistic to collect statistical data on biomass for good size data sets. 

It is unrealistic to expect statistical data on the number of leaves and 
their orientation, twig volume, and number of flowers or fruits for anything 
but a small collection of cases . Due to the scarcity of good statistical 
data on target characteristics (i.e., biomass, water content) reflectance 
variability was driven by an artifical variable. Statistical data at the 
process resolution desired was not always available or complete and estimates 
were used where necessary. For targets, statistical data was not available 
for the fundamental characteristics or for the reflectances. 

6.2 Statistical Distribution Characteristics 

The assumption of the normality of statistical distriubtions in remote 
sensing work is almost universal. For those quantities, such as sensor 
irradiance, which are a function of many independently (perhaps) distributed 
quantities, the central limit theory of statistics would provide some jus- 
tification for that assumption.. 
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Strictly speaking, the normal distribution is not valid for quantities 
such as path radiance because of the exclusion of . negative ■ radiances . The 
normal dis.triubtion was chosen for the distribution of the atmospheric con- 
stituent amounts and reflectance artificial variables. For atmospheric 
constituent amounts this assumption is intuitively pleasing. The case of 
reflectance variables does not offer such appeal. One problem concerns 
the rather arbitrary nature by which targets are identified. For instance, 
one may ask at what point does thin vegetation become bare land or what 
is the difference between a thin cloud and thick haze. If for example, the 
reflectance variable for vegetation were biomass, there seems no intuitive 
justification for assuming that it is normally distributed. A more likely 
case would be an asymmetrical distribution that peaked at the biomass deemed 
to be the necessary minimum to be classified vegetation. The second major 
limitation with regards to reflectance artificial variables is the restric- 
tion of reflectance to the range of 0 to 1. For a particular portion of 
the spectral reflectance curve for a target characterized by low or high 
mean reflectance, the range of possible reflectances is truncated. If the 
normal assumption is used for the artificial variable in this situation, 
the allowable size of the reflectance variability is unduly restricted. 

The failure of the normal assumption for radiance data has been 
demonstarted both in field and computational experiments. Valovcin (1978) 
plotted mean radiance and ±1 sigma radiance curves for snow and cirrus 
clouds. Both of the minus sigma curves showed radiances below zero which 
are physically impossible. Presenting ± sigma curves (or bars) suggests 
things about the statistical distribution that may or may not be true. 

Sigma curves, which are symmetric about the mean, lead the reader to believe 
that the probability density function (PDF) is symmetrical. The plots 
presented by Valovcin clearly preclude the assumption of a symmetrical PDF 
and therefore rule out the normal assumption. Similarily the 1 sigma 
ellipses for several simulated classes in Fig. 4.1 enclose negative radiances 
in the same manner as Valovcin' s data. 

6 . 3 Distribution Parameters 

It was mentioned briefly in the section on reflectance there is a cor-, 
respondence on covariance between channels and their redundancy. The exact 
redundancy is a function of the probability density functions of the channels 
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For a multivariate normal distribution the channel redundancy is speci-- 
fied exactly by the covariance matrix (see for example a text on Infor- 
mation Theory) and this is the distribution that will be assumed for 
illustrative purposes. Using the same normality assumption, Wiersma and 
Landgrebe (1979) attributed all the covariance for a set of spectral radi- 
ance data to covarying reflectances, with one exception. The atmospheric 
water vapor absorption bands were eliminated from the analysis. The 
sophistication and rigor of their analysis is almost unparalled in the 
remote sensing literature, and should be reviewed to understand the impor- 
tance of covariance "structure" in determining optimal spectral channels. 

One of the shortcomings of Wiersma and Landgrebe ’s analysis is the lack of 
consideration given to the contribution of the atmospheric variability to 
the spectral radiance covariance structure. This shortcoming is almost 
certainly due to the lack of data to quantify and statistically characterize 
the atmospheric effects present in the data. 

The simulation that has been discussed herein has the capability to 
combine the effects of the reflectance and the atmosphere on the spectral 
radiance covariance structure. The reflectance covariance structure was 
discussed in the surface reflectance section. Defining the atmospheric 
contribution to radiance covariance has been no more fruitful than for the 
reflectance. The model utilized for atmospheric radiative transfer vari- 
ability requires that the covariance of atmospheric constituent amounts be 
defined. For example, it would be necessary to define the correlation 
between atmospheric water vapor content and aerosol burden. So far no 
such data has been found. Simulations run to date have assumed all consti- 
tuents to be uncorrelated with the exception of aerosols and water vapor 
which are assumed to have perfect positive correlation. Fraser (1975) re- 
ports covariance data for optical depths determined from solar attenuation 
data. By assuming some model for the various atmospheric constituents 
contributions to optical depth, it might be possible to invert the covariance 
data for the constituent amounts but so far efforts in this area have not 
been successful. 

6.4 Decision Theory Considerations 

It became apparent during the work represented in Aherron, et al. (1981) 
and cloud detection algorithm performance reported herein that classification 
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or categorization accuracies where highly dependent. on the probability of 
encountering various targets. This is a familiar subject to those with 
knowledge of Bayesian decision theory. -^The two statistical values necessary 
for Bayesian decision calculations are the conditional probability of a 
signal given that it was produced by a certain class and the a priori pro- 
bability of encountering that class (Duda and Hart, 1973). The statistical 
characterizations and models discussed to date have been concerned with 
defining the conditional probability function for a class. The algorithm 
performance evaluation model utilized thus far has been a simple overall 
classification or categorization accuracy. These figures have been based 
on equal a priori probabilities for classes. When classes are aggregated 
into categories the a priori probabilities are not necessary equal. This 
point must be kept in mind when comparing the reported overall decision 
accuracies. 

Discussion of the Boundary Approximation Method (BAM) decision 
algorithm in Aherron, et al. (1981) points out the need for class a priori 
probabilities in balancing the type I and type II errors of categorization 
when formulating the exact algorithm parameters. As the overall remote 
sensing model becomes more realistic and accurate in future work, it will 
be necessary to include realistic a priori probabilities of occurence for 
both classes and categories. At least two potential methods for channel 
bandwidth and location selection will require such knowledge. One of these, 
an information theory based method, was used by Kondratyev (1975) for 
choosing spectral channels. Kondratyev's analysis was fairly crude and 
the method he used is still relatively unexplored. Secondly, Bayesian 
Decision Theory, as was discussed before, requires the same class 
probability data. 
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7. SUMMARY & CONCLUDING REMARKS 

This report has presented supporting research and data for a previously 
reported comprehensive remote sensing system simulation (Huck, et al. 1982). 

A review of various models of atmospheric radiative transfer was presented. 

The qualitative accuracy and the quantitiye computational complexity of the 
various models was compared. Using those comparisons and evaluation criteria 
that included recognition of the uncertainties about atmospheric parameters, 
the Turner (1974) model was chosen for use in the remote sensing system simu- 
lations. A detailed discussion of the stochastic implementation of the Turner 
model was given. Various functional forms for spectral diffuse reflectance 
that allow introduction of variability were given and discussion was made of 
certain properties of reflectance variability that are desirable in dealing 
with realistic targets. A compilation of various spectral reflectance curves 
utilized in simulations was presented. 

The special problem of detecting clouds reliably has been examined. For 
optically thick clouds a test of a preliminary detection algorithm has shown 
very high accuracy. Thin clouds pose a problem in detection that needs to be 
examined more closely. The detection algorithm tested requires a spectral 
channel at .65 ym and 1.55 pm. As was pointed out in Huck, et al., (1982), 
the remote sensing system simulation being developed is unique in its treatment 
of the stochastic elements of the remote sensing system. The lack of good data 
for parameterizing the statistical nature of the atmosphere and target reflect- 
ances was discussed especially as it pertained to defining the wavelength covariance 
matrix for radiances. The need for a priori probability of occurence data for 
targets was discussed. Finally, a partially updated set of categorization 
results for the experiment reported in Aherron, et al. (1981) was presented 
which showed the simplest method of categorization performing even better 
relative to the other methods than previously reported. 

Hopefully modeling efforts such as described proyide insight to the over- 
all sensor systems design effort as well as those people involved in basic 
research concerning fundamental processes. 

One role that remote sensing system modeling plays in research is to 
encourage a more Insightful look at ayailable data and more careful defini- 
tion of needs for future data collection. Specifically the statistical nature 
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of the various system elements becomes, more important as the modeling 
becomes more refined and generalized. And finally, modeling should help to 
make informed decisions as to the importance of various effects and error 
sources within remote sensing systems. 
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